Noise induced stabilization in population dynamics 
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We investigate a model where strong noise in a sub-population creates a metastable state in an 
otherwise unstable two-population system. The induced metastable state is vortex-like, and its 
persistence time grows exponentially with the noise strength. A variety of distinct scaling relations 
are observed depending on the relative strength of the sub-population noises. 



The phenomenon of noise induced metastability 
0| is of importance in ecology [1] and plant biology _ 
and has found practical applications in engineering 
The typical models f]]-!^ consider periodically modulated 
one-dimensional (Id) stochastic systems. The modula- 
tion renders the system to be deterministically unstable 
during a part of the modulation period. An external 
noise can prevent the escape for several successive peri- 
ods of external modulation, trapping the system into a 
metastable state. As a result, the noise causes an increase 
of the system persistence time by a factor compared to 
the noiseless case. 

Here we consider a different model with two stochastic 
degrees of freedom, which we call x and y. The y de- 
gree of freedom (e.g. imbalance between the numbers of 
two competing gene alleles) undergoes a strong and fast 
noise which conserves the total population size x. The 
latter experiences a slow evolution under the influence of 
a deterministic potential V(x) along with a sign-definite 
feedback from the population size imbalance oc y^ and 
a relatively weak (demographic) noise. We show that, 
even if the x-dynamics itself is unstable and prone to 
a rapid escape, the strong y-noise can lock it in an ex- 
ponentially long-lived vortex-like metastable state. The 
corresponding exponent exhibits a variety of non-trivial 
scaling regimes, depending on the relative strength of the 
noises in the x and y subsystems. A similar model was 
shown to describe a two-patch Lotka-Volterra system . 
More distantly related models were recently discussed in 
the context of biochemical regulatory networks [s*] and 
nanomechanical oscillators 

Our model can be cast into the universal form 



x = -V'ix)-y^+Ut), 
il^-2y + iy{t)- 

{ix(y){t)ix(y){t')) = '2T^{y)5{t - t') , 



(1) 



where V — dV{x)/dx, and and Ty characterize the 
noise strength in the total and differential population 
size, respectively. The interesting regime of parameters 



is 



< Ty. The noise effects are substantial when the 



population is close to a bifurcation point. In this case the 
properly rescaled deterministic potential takes the form 



>> 1.0 




FIG. 1. (Color online) The quasi-stationary FP state. The 
contours represent the probability density P{x,y); the ar- 
rows show the probability current density. The system is 
symmetric around y = Q, and only the top half is shown. 



= 0.05, T„ 



0.5. 



V{x) = -x^/i - 5x, 



(2) 



where 5 is the bifurcation parameter, and we have shifted 
the X variable to have the bifurcation point at a; = 0. 

A simple realization of this model is provided by two 
species A and B, which undergo the reaction A-\- B \ 
2A, where X is either A or B. This is the well-known 
Moran process for modeling neutral genetic drift [lo|. 
This is the fastest process which conserves the total 
population size. In addition, the population size may 
slowly evolve according to e.g. the following set of reac- 
tions A and A + B A + B + X. In this case 
X = [ua + — N)/N and y — (n^ — nB)/N, where 
N ~ f3^ / a is the population size close to the bifurcation, 
and 6 — 4a/3+//3^ — 1 is the bifurcation parameter. In 
the limit of large population size N ^ 1, the correspond- 
ing Master equation can be approximated by a Fokker- 
Planck equation ll|. In the vicinity of the bifurcation 
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point, the latter reads 



P{x, y) =-d^ 

-a 



{-V'{x)-y^) P{x,y)--d^P{x,y) 



-2yP{x,y)- [\+-]dyP{x,y) 



(3) 



where P{x,y,t) is the probabihty distribution function, 
and time is measured in units of 2/I3-. Equation (|3]) 
is equivalent to the Langevin equations ([ij, where the 
two "temperatures" are given by = 2/N and Ty = 
Tx + 2A//3_. When the drift rate A is fast, one has the 
strong inequality Ty. 

First we focus on the case of exact bifurcation, 6 — 0. 
The x-equation takes the form x — x'^ — y"^ + With- 
out noise the y-variable tends to zero, leading to i = 
dynamics in the x-direction. This has a; = as the 
marginally stable point. An arbitrarily weak a;- noise is 
sufficient to kick the system out of this fixed point and 
set it on the path to unlimited proliferation, x — >■ oo. 
One may think thus that the ^ = system is destined to 
blow up in a very short time. Recall, however, that the y- 
noise is substantial. Although (y) = 0, the mean square 
value (y^) > and is large compared with T^. One can 
then expect the x-dynamics to be governed by the effec- 
tive potential Vcff(x) = — x'^/3 -I- (y^)x. This potential 
exhibits a minimum at x = — y'Jy^, and a maximum at 
X — \/ (y^). As a result, a long-lived metastable distribu- 
tion, peaked at x = —\J (y^), can be created. A numeri- 
cal solution of the Fokker-Planck (FP) equation ([3|) sup- 
ports this expectation. Figure [T] shows the slowly vary- 
ing quasi-stationary distribution observed at late times. 
Notably, the probability currents develop two counter- 
rotating vortices. Before reaching the point x = y = 0, 
the "particle" is kicked in the y-direction, where the x- 
evolution is directed toward population contraction. The 
vortices, therefore, arrest the population explosion. We 
note that probability current vortices in non-equilibrium 
stationary states ~ the Brownian vortices - were recently 
observed in experiment 12]. 

Our main goal is to evaluate the lifetime of such a 
noise-induced vortex-like metastable state. We start 
from qualitative considerations. As a first approximation 
one can estimate the mean-square y-deviation in the har- 
monic potential y2, cf. Eq. (P), as (y^) = Tyjl. Therefore 



T/max 



''off 



\f2Ty /3, and one expects that 

iniesc=^ V2r3/V3r,. 



(4) 



Remarkably, the escape time is exponentially increasing 
with the y- noise strength Ty, while exhibiting the stan- 
dard Arrhenius scaling with T^. Our numerical simula- 
tions of the Langevin Eqs. ([T]), see Fig. [21 confirm Eq. 
as long as Ty^^ /T^ is not too large. At larger values of 
this parameter, however, Eq. ^ greatly overestimates 
the lifetime of the metastable state. 
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FIG. 2. Simulated escape times for Ty = 0.05 and varying T^. 
The straight line has slope of \/2/3, cf. Eq. Q. Inset: the 
limit yJTZ > Ty. 



The reason for this deviation is that for large Ty the 
typical potential barrier is too high for the x-motion to 
overcome. Then, instead of relying on typical realiza- 
tions of y-noise, the system prefers to wait for a rare y- 
trajectory which stays anomalously close to y = 0. The 
probability that the y-motion is confined to the inter- 
val \y(t)\ < yo for a time to is given by exp[— ii^y (yo)to]- 
Here Ey is the lowest eigenvalue of the Id FP equation 
in the y-direction with absorbing boundary conditions at 
y = ±yo. Ey can be estimated as Ey cx Ty/y"^, where 
Ty is the proper diffusion coefficient. On the other hand, 
the probability that during the time interval the x- 
coordinate will diffuse from x = — yo to x = +yo is given 
by exp(— yQ/Tj-tg), where T^ is the proper diffusion coef- 
ficient. Maximizing the product of these two probabili- 
ties with respect to yg/toi one finds that the probability 
of the optimal rare fluctuation scales as exp{— ^^Ty/Tx). 
These estimates suggest that Iniosc oc ^TyjT^ once 
^TyjTx < Ty^^/Tx, i.e. -v/TT < Ty. This behavior is 



indeed qualitatively consistent with Fig. [5] 

To put these considerations on a more quantitative 
basis we shall assume that the dynamics can be sepa- 
rated into the fast y-motion and slow x-motion. The 
latter adiabatically adjusts to the instantaneous value 
of y^(i). We then solve an auxiliary problem of find- 
ing the probability of y-trajectories with a given func- 
tional form (y^) — y^it). Here yo(t) is an arbitrary 
slow function of time, such that yo(±oo) = ^jTy/2, 
while the averaging is taken over the fast y-fluctuations. 
Integrating over an intermediate time-scale At that is 
long relative to these fluctuations one can then write 
It-Atlyoi^) ~ y'^{'^)]d't = 0. We can thus introduce 
the functional constraint i5 ( / [y^ {t) — yn (tlldt) into the 
stochastic functional integral over Vy |13l4l5l ] and elevate 
it into the exponent with the help of the auxiliary slow 
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FIG. 3. Phase portrait of the optimal escape paths for differ- 
ent Tx. Lower lines correspond to lower values of T^. The 
momentum has been rescaled by Ty/Tx so that all paths 
would coincide if they followed simple activation [assumed 
by Eq. @]. 



FIG. 4. Simulated escape times vs. the optimal action for 
Ty = 0.05 and varying T^. Inset: the limit {T^/T^f^'^ > 1 
(with Ty = 1000). The straight line is Eq. 0) . 



field x(t). As a result we obtain an effective Lagrangian 



r — 

l^y — 



4r„ 



X(y^ - Vo) . 



(5) 



where the ^-integration runs from —ioo to ioo. Employ- 
ing the slowness of the x(i) field, the Gaussian integral 
over the fast y{t) can be evaluated using the Fourier 
transformation. This leads to an effective Lagrangian 
for X in the form 



xvl 



dco 
2^ 



In 1 



^TyX 



xvl 



1 + 1-^%^ 
(6) 

Finally, the x-integration can be evaluated in the saddle 
point approximation: x{t) = Ty/Ay^ — Ty^. This yields 
the probability of y-niotion conditioned on (y^) = y1{t): 



PMoce-^*^'''^^), Ey{yo) = ^ 



l + (7) 



Notice that Ey is non-negative and equal to zero if and 
only if j/q — Ty/2. Therefore, the condition ?/g(±oo) — 
Ty/2 is necessary for convergence of the integral in 
Eq. ([7]). The saddle point calculation is justified as long 
as JdtEylyoit)] > 1. 

Having found the conditional probability of y-motion 
with a given profile of (y^), we turn now to the x- 
degree of freedom. According to the scale separation 
assumption, it is governed by the Langevin equation 
i — — y^it) + £,x{t), where yQ{t) is a slow function 
of time with ?/g(±oo) = Ty/2 and yg(0) < Ty/2. Our 
goal is to evaluate the escape rate of the x- variable from 
its metastable minimum at x — —^Ty/2 during the 
time when y\ it) is suppressed with respect to its asymp- 
totic values. We then maximize this escape rate, taken 
with weight P[?/o], Eq- 0, against the optimal time- 
dependent variance yo(i)- 

Since the escape rate in the x-direction is expected to 
be small, it can be found through a semiclassical treat- 



ment of the corresponding FP equation 16| . The proper 
FP Hamiltonian has the form 

U\x,vx-y^{t)\ ^p^\-V\x)-yl+TxV.\~Ey\y^{t)\, (8) 

where x and are canonically conjugate variables, and 
2/0 (i) is an external time-dependent parameter. The last 
term accounts for the statistical weight of a realization of 
2/0 (i), given by P[2/o]j Eq. ([7]). If yo(^) is an adiabatically 
slow function, the escape proceeds along the zero-energy 
trajectory of this Hamiltonian, which connects the two 
fixed points (-^7^2,0) and (+\/Tj^/2,0) on its {x,Vx) 
phase plane, Fig.|31 Putting V{x^ — —x^/Z, one finds for 
the (slowly varying in time) optimal trajectory 



1 

2Tx 



+ ^{yl-x^Y+ATxEy{y^) 



Px{x\yo) 

(9) 

The corresponding escape time, within exponential accu- 
racy, is given by the classical action, i.e. the area of the 
phase plane under the zero-energy trajectory 



lnfesc[j/o] = 5' [2/0] 




Px{x;yo)dx. 



(10) 



The final step is to find the optimal 2/0 realization. This 
is achieved by demanding 5S[yi)]/5yQ = 0, solving for an 
implicit function of time 2/0 = 2/0(2;) and substituting it 
back into Eq. (flUl) . This leads to the optimal action, 5opt , 
and corresponding escape time Iniosc = 'S'opt- In Fig. |4] 
this escape time is compared with our Monte-Carlo sim- 
ulations results, and an excellent agreement is observed. 



It is easy to show that, for Ty <^ 



2/0 tends to yjTy/2 and thus Ey 



V2Ty'^/iTx. 



We thus recover Eq 
site limit Ty » one finds 2/0 (0) < \/Ty/2 



/Tx, the optimal 
0, while ^opt = 
(gl). I n the oppo- 
One 



can thus simplify Eq. ^ as Ey k, Ty/Ay^. With this 
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FIG. 5. Phase diagram of the system as a function of and 5 
with Ty held constant. The dashed lines represent crossovers 
between different scaling relations. In region I, Infcsc = 
In region II, Infesc = 4 {Ty/2 - Sf'^ /ST,. In 
region III, In icsc = 7rTy/2(5^/^ 



substitution Eq. can be simplified by the rescaling 
X = iiT^Tyf/^ and yo = yoiT^Ty)^/^, which brings the 
action (llOp into the form S = y^Ty/T^ J p(x; yo)dx. The 
integration limits are ±{T^ /TxY^^ ±oo in the limit 

' „--.-2 



of interest, while p ■ 



/2 



is a parameterless function. Optimizing it over yo and 
performing ai-integration, one finds 



In 



^T, « Ty 



(11) 



which confirms our qualitative estimates below Eq. ^ 
and provides the numerical factor. The latter is com- 
pared with our Monte-Carlo simulations in the inset of 
Fig. m Notice that the actual condition for the applica- 
bility of the asymptotic result ^ is 1 < {T^/T^)'^/^. 
Again, the population lifetime increases with the y-noise 
strength. Notice also that the normal Arrhenius scaling 
in the parameter Tx gives way to a stretched exponential 
law with Tx~^^'^ . A similar transition is known as Efros- 
Shklovskii law [r?| in the context of hopping transport in 
disordered semiconductors. 

We consider now deviations from the exact bifurcation 
point, i.e. 5 ^Q. If \5\ ^ Ty, the deterministic "force", 
cf. Eq. is very strong, and the y-noise is not impor- 
tant for the system's persistence time. Pd 5 = 5c — Ty/2 
the effective force associated with the y-noise is canceled 
by the deterministic i5-force. This causes a noise-induced 
shift in the bifurcation of the cc-dynamics. That is, it 
is much harder to destabilize the population in the pres- 
ence of strong y noise. In the vicinity of this noise-shifted 
bifurcation one finds the standard scaling of the lifetime 
In W = 4 ((5c - Sf'^/^Tx, cf. Eq Q. 

On the other hand, away from the the noise-shifted 
bifurcation, i.e. at \5 — 5c\/&c > iVT^/Ty)'^^^, the scaling 
changes qualitatively. To find the new scaling we look for 
the zero energy trajectory of the Hamiltonian (|S]) with 



6 and Ey = Ty/Ay^ and optimize the action over 
yo(x), as explained above. In this way we find 



In U 




[TxTyY/^ 



'Tx<^Ty, (12) 



where the universal function S{5) has the following 
asymptotic limits: S{5) w 1 -I- 0.715 for (5 <C 1 and 
S{5) « 3 (5^^/'^/4 for (5 ^ 1. This means that the scaling 
of the population lifetime given by Eq. ([Tl]) is basically 
intact as long as 5 ^ {TxTyY^^. In the opposite limit, 
the escape time scales as 



Intcsc = ■nTy/25 



3/2 



(13) 



This is independent of T^. For 5 > [TxTy^/^ the system 
can escape even at T^; = via paths with unusually small 
y. Figure [5] shows the regions where the three scaling 
relations g]), ([11]), (HH) are vahd. 

In summary, we have studied a novel system where 
strong noise creates metastability. Increasing the noise 
strength Ty increases the lifetime of the (vortex-like) 
metastable state. Escape from this state is governed by a 
variety of scaling relations depending on the relative role 
of Tx (the strength of noise in the total population size) 
and Ty (the strength of noise in the differential size). 
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